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ABSTRACT 

In this work, variations on velocity profiles in a flowing mass 
of molten glass in a forehearth are investigated. Formerly used 
parabolic velocity profiles are replaced with analytical solution of 
open channel flow equation, based on the available data on mass flow 
of molten glass through the channel in unit time. 

Concerning the viscosity effects; temperature dependence of 
viscosity is built in the model. However, it is assumed that the depth 
of the channel does not affect the viscosity gradients. 

To solve the system of non-linear differential equations i.e., heat 
equation and flow equation; the analytical solution of the latter at 
the nodes is used for the numeric solution of the former iteratively, 
until the convergence is obtained. Predicted temperatures are compared 
to the available data from an actual operating forehearth, and against 
the results predicted by the previous model using simplying assumptions 
to prove its validity. 
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I. INTRODUCTION 



To improve the production of glass containers, and to transform the 
art of glass making into a science, the glass industry embarked an 
extensive forehearth measurement program in an attempt to resolve the 
complex phenomena of glass conditioning and cooling, while molten glass 
is taken out of the furnace and transferred to the forming machine. 

The objective of the program was the collection of quantitative data 
to achieve greater understanding of forehearth and to verify the math- 
ematical model of molten glass flow in forehearth developed by Duff in 
[Ref. 1]. 

A knowledge of the temperature distribution is needed for better under- 
standing of the control process in glass plants. The above model proved 
that sophisticated models can give satisfactory results describing the 
forehearth behavior in a temperature sense. However they need a big 
computer and require considerable time to run. 

In this work, validity of the simplifying assumptions in [Ref. 1] 
are investigated and ways to relieve them are sought. Going to a more 
complicated model has academic interest to show that what was developed 
was a useful and fairly accurate model. 
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II. NATURE OF THE PROBLEM 



Forehearth is that section of the process where molten glass is 
transferred from furnace to the forming machine. During the flow, glass 
is conditioned in a temperature sense. A predetermined temperature is 
desired at the forming machine end, which assures minimum off-grade 
material after forming process. If undesirable temperature gradients 
exist at the input of forming machine, recycling of off-grade glassware 
is necessary, decreasing the overall system efficiency. 

What is needed is controlled conditioning throughout the forehearth 
while molten glass is flowing. This is mainly a heat transfer problem 
on a flowing mass of molten glass. 

To describe the physical phenomenon taking place in a typical fore- 
-hearth used in glass container manufacture, see Figure (1). Glass flows 
from left to right in a open channel with dimensions approximately 26 
inches wide, 6 inches deep and 18 feet long. 




Section AA Section BB Section CC 
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Section M is the point of entry from furnace area with average 
temperature being 2400 degrees F. at this point. Area from AA to BB 
is called cooling zone. Area from to ^ is known as conditioning 
zone. This is the place where an attempt is made to condition the glass 
so that the delivery temperature would be constant in plane After 
CC a glass gob is formed. A gob is a discrete mass of molten glass 
created by intermittenly shearing a stream of molten glass coming 
from an orifice. The gob formed falls into a guide chute, and is con- 
ducted to a forming machine mold. 

Temperature control is achieved by adjusting burner flame-levels and 
the amount of cooling wind in the cooling zone (section M to ^); also 
with burner flame levels in the conditioning zone. 

Burners are mounted on manifolds and are located every 4.5 inches 
on both sides of the forehearth. They produce a blanket of heat over 
glass surface. Hot gases sweep over the ceramic radiating surface and 
heat exchange takes place between gases and surface which, in turn 
exchanges heat with molten glass. Cooling is achieved by introducing 
cooling air through inlet pipes. The wind exchanges heat with glass and 
the ceramic surface and goes out of the system through roof ports. Details 
of the forehearth cross section is shown in Figure (2). 

Heat exchange at the ceramic surface takes place by convection and 
radiation mechanisms. Heat is also conducted through the ceramic material 
in the sides and bottom of the forehearth. Possible disturbances can 
affect the system through these walls. 

Above discussion indicates a complex problem of heat exchange in fore- 
hearths, with boundaries as described. 

To describe temperature behavior within the glass, it is necessary to 
consider the mechanisms by which heat is transferred. 
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It had been shown that emission and absorption of radiation is a 
bulk phenomena. The interaction of simultaneous emission and absorption 
throughout the volume leads to radiative heat transfer between adjacent 
layers of material. When this radiative heat transfer is combined with 
true conductivity, mathematical description of this internal radiative 
mechanism becomes rather involved. 

Glass is, then, a diathermanous material whose diathermancy is 
determined by 

a) Wave length-absorption coefficient relation and 

b) Spectral composition of radiation involved. 

This dependency is usually shown in plots of absorption coefficient y 
versus wave length A, with temperature as a parameter. Any detailed 
analysis of the internal radiation should recognize this varying 
diathermancy. Integration can be replaced by summation using average 
values of y for various ranges of a, with the average value in a given 
A range being a function of temperature. Because of its transparency, 
glass layers beneath the surface will exchange heat with the surroundings. 
The depth to which this exchange occurs is appreciable varies with 
the type of the glass being considered. Upper one third of the glass flow- 
ing in a forehearth is considered as exchanging energy directly with the 
surroundings by Duffin [Ref. 1] in his model. 

In the glass body, the radiant effects result from bulk phenomena of 

emission and absorption. Internally emitted radiation is reabsorbed by the 

glass directly and also as a result of internal reflections. This process 

leads to the "radiative conductance" of heat through the glass. Internal 

emission is characterized by '.'volume emissive power" which is the rate at 

which a unit volume of glass emits radiation in all directions. For an 

ideal gray material volume emissive power is given as, 

W = Ayn^oT^ 
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ThuSj it is proportional to absorption coefficient, refractive 
index and Stefan-Boltzman constant. According to Kellet [Ref. 2], in 
the interior of massive bodies of molten glass, steady-state temperature 
distributions turn out to be linear. It was also proposed that the 
effect of radiative conductivity be designated as an "radiation conduc- 
tivity". For the gray material it is given by 

_ 16 n^aT^ 

•^RAD “ 3y 

Thus, the effective conductivity within the material is the sum of the 
thermal conductivity and radiation conductivity, as follows: 

^ " *^eff " “^true * “^rad 

This is set and held throughout the glass depth. But due to 
dependence on the absorption coefficient, variations are possible near 
the glass surface. 

A. DERIVATION OF DIFFERENTIAL EQUATIONS 

Derivation of the differential equations was accomplished by Duff in 
iRef. 1] as follows; using a right handed coordinate system and taking 
X coordinate as the flow direction, a differential volume element of 
dimensions dx, dy, dz is set up. Then writing energy balance with 
dz = 1. 

1 . Energy Into Element Due To Mass Flow 
Area in X direction = dydz = dy 
Mass flow = V (1 *dy)p 

Energy flow = (pdyV^)(C)(T) 

2. Energy Out Of Element Due To Mass Flow 

| 3 ^ (pCVJdy)dx + pCVJdy 
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3. 



Net Flow Energy 



= Input - Output = 

- ~ (pCV^Tdx)dy 

4. Energy Into Element Due To Conduction-Radiation 

Area = dx-dz = dx 
gT 

Energy = - k'dx — 

oy 

5. Energy Out Of Element Due To Conduction-Radiation 

|y{-k'dx|I)dy.[-k'dx|l] 

6. Net Conduction-Radiation Energy 

I7 

7. Net Total Energy Flow 

ly I?) - k (pCV,T)dxdy = (dxdyp)C- |I 
sOs differential equation is 

This being one dimensional model, actual case must include the wall 
effects which is not considered in above derivation. To account for that, 
a term for Z direction must be included. 

I7 ly) ^ I? It) - k (V'T) = pc |i 

gT 

If — is equal to zero, steady-state heat flow equation is obtained. 

In above equation density, specific heat and effective conductivity terms 
can be treated as temperature independent. Dependence of velocity on 
temperature will be discussed later. 
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B. BOUNDARY CONDITIONS 

To apply boundary conditions, a fixed coordinate system is placed 
in channel center as follows 




T 


(3 


X = 0 


y and z = 0 


Temperature in the initial plane 


T 


0 


y = d 


X and z = 0 


Radiating boundary 


T 


0 


y = 0 


X and z = 0 


Bottom temperature distribution 


t 


0 


z = w 


X and y = 0 


Side wall temperature distribution 



at y = 0 and z = w glass contacts with the ceramic material. The tem- 
perature of the glass and ceramic material will be same at these points. 
Following will also apply 



(il) 

'3yV=0 



c 

F" 



(— ) 
^ 3y^ 



=- _£ I 9.) 

'32^Z=W k' ' 3Z' 



The boundary at y = d is the radiating boundary. In the model it is 
assumed that the enclosing surface and the glass surface behave as two 
infinite opposed parallel planes. If h is the equivalentx^coefficient of 
heat transfer for the convection-radiation exchange for the volume element 
below 




A 

conduction 
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Energy Input 



Area = dx'l 



#y=6 * 

Energy Output 

hdx (T - Tg) + aegT^ dx 
Accumulation 



-k'dx dx-hdx(T - T^) 



ae«T dx = pC^6dx 

6 p 
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Let 6 -»• 0 , so y -»■ d 

#)y=d ■ h f - t') - ^ (T - T^) 

is the glass to air enclosure boundary condition. F is the view factor 
which accounts for the geometric arrangement of the two surfaces. With 
the assumption above following equation applies: 

F = ^ 

l/e^ + 1/eg -1 

Where emissivities of the enclosing surface and 

glass, respectively. To go one step ahead, this definition of view 
factor is replaced with a new one, which also accounts for the re- 
radiation from the connecting walls. Glass and ceramic surface being 
gray surfaces connected by nonconducting, reradiating walls, the new 
shape factor is given by Chapman [Ref. 3] as 



F = 



1 



+ ^ (i_ 
fl-2 *2 'S 



- 1 )+ 1 ) 



A - A F 

^2 n n-2 



1-2 






2 A.|F^_2 
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F-j 2 IS the geometric shape factor, A-j and A 2 are the areas of the 
radiating surfaces. 



C. NUMERICAL METHOD USED TO SOLVE THE HEAT EQUATION 

Partial differential equation describing the steady-state temperature 
T(x,y,z) in a moving rectangular slab of molten glass in a forehearth is 



+!<1. A-.9- (VT) = 0 



( 2 . 1 ) 



Where cartesian coordinates are used,x is the flow direction down 
the length of the slab from the front of the hearth, y is the direction 
from the bottom channel toward the top of the slab, and z is the direction 
from the channel center toward the side, k', k", Cp are the known 
constants. The initial temperature distribution at x = 0 is assumed 
to be linear in y and z. It can be found by linear interpolation from 
the four corner temperatures. Similarly it is assumed that the temper- 
ature on the boundaries (y = 0), bottom of the channel, and (z = w), 
side of the channel, are also linear. I.e., 



T(0,y,z) = <|>i(y,z) 


(2.2) 


T(x,0,z) = (|> 2 {x,z) 


(2.3) 


T(x,y,w) = c|) 3 (x,y) 


(2.4) 



and are all known. The temperature of the glass is symmetric with respect 
to the center line (z = 0); therefore. 

Radiative boundary at the surface of the glass (y = d) is that of 
radiative and convective heat transfer. Equation derived previously was. 
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Using linearizing approximation on the nonlinear part, and changing 
sign on the second term of right hand side, 

(|y)y=d ■ ¥ - Vd> * ^ - Vd> <2-6) 

Where o, F, k', h are known constants. T, is the air temperature on the 

a 

glass surface. T^, the temperature of surroundings is determined by 
linear interpolation between the temperature at the center of the channel, 
T$o and the known temperature on the boundary at the side (z = w). The 
center temperature is assumed to decrease linearly from initial temper- 
ature T at X = 0 by an amount T , at the end of the channel (x = XL); 

Si c 

therefore 

= ■'si - IT T, 

With these, the problem is formulated as being the solution of the 
partial differential equation (1), with initial condition (2.2), and 
boundary conditions (2.3), (2.4), (2.5), and (2.6). 

An unconditional ly-stable alternating-direction method [Ref. 4] was 
applied after passing equation (1) into finite difference form. This 
is accomplished by first getting the first derivatives with backward 
differences then substituting these in forward difference equation to 
get finite difference equation for second space derivatives, in y and 
^ directions. Forward difference form is used for the derivative in x 
direction. 

The first increment step is taken as implicit in y direction, the 

second implicit in z direction, and so on, x step size being the same in 

each case. Temperatures at successive planes in x direction are related 

to each other since the first step involves values at x + ax and these 

*n+l 

intermediate values (T. . )are used in the following step during the 
numeric calculation. System of equations obtained are. 
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^i.J> - '^l<^i-l.J - 2Ti,j T,+I,j) + C2 (T,_j.,-2T,^j+T,_j^,) 



r*n+l 



*n+l 



- T*"j^= C,(T*n+l _ 2 T*n+l + -r*n+l x. . (Jn+^ .,n+l, .n+1 ^ 

With following definitions using increments ax, Ay, Az 

bij = V^(iAx,JAz) 
r _ 1/2 ax k' 

’’IT" ^ 




k" 




t!? j = T(nAx, iAy, Jaz) 

if 

Ti,j = T(nAx, iAy, Jaz) 

The indice i runs from i = 1 to i = I where i • Ay = d, indice j runs 
from j = 0 to j = J-1 where j*Az = w. The initial temperature gives 
the values tV ^ = 4 i.j(y,z), boundaries give the values Tq j = Tq j = 
(f> 2 (nAx, Jaz) when i = I and t!? ,j = T**^j = c)) 2 (nAX, iAy) when J = J-1 . 

Also, symmetry around the center line of the channel gives rise to 
t|? 1 = T^ , and T.*^ ^ = T. , when j = 0. Radiating boundary results 

in substitutions for t!?^i j and j in terms of T^ j and T^*^j when 
i = !• 4>i . 4>2 4 >^ are obtained from corner temperatures by linear 

interpolation. 

The system of equations is solved starting from t 9 ., getting the 

T sJ 

intermediate solution from the first and using it in the second, making 
use of Thomas algorithm [Ref. 5] to solve resulting tridi agonal matrices, 
and proceeding through the channel at every y, z plane separated by 
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AX in X direction, until the end of the test section is reached (x = XL). 
For integration purposes x length of the channel is divided into 500 
increments (ax = XL/500), y coordinate is divided into 30 increments, 
and z coordinate is divided into 75 increments, making use of the sym- 
metry around the center line. 
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III. VELOCITY PROFILES 



A. PARABOLIC 

First considerations of velocity distribution within the flowing 
mass of molten glass by Duffin [Ref. 1] had involved moving side walls 
of the channel to infinity which gave velocity in x direction as a 
simple function of y coordinate. Next step was evaluation of finite 
velocity step surfaces based on a given mass flow rate of glass in an 
effort to account for higher velocities in the center channel and lower 
velocities near the side walls and bottom. None proved satisfactory , 
so final choice was the use of explicit expressions for parabolic 
velocity profiles. Physical picture is as follows: 




Maximum velocity occurs at y = d and on center line z = 0. Derivation 
of the general equation of parabola families gives velocity in x direction 
as a function of y and z, and physical parameters w, d and as follows 



V (y,z) = V 
x^-^’ max 






( 1)2 + (^)2 
vw^ ^ ^dw’ 



(3.1) 



from which 9/4 = V„,^ can be obtained. After integration over the 

qV6« iTIqX • 

channel cross section and dividing by channel area. 



B. VELOCITY PROFILES USING OPEN CHANNEL FLOW EQUATION 

Above parabolic profile does not account for the viscosity-temperature 
relationship. Cooper [Ref. 6] pointed out that the effects of moderate 
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viscosity gradients on velocity are negligible. Temperature in the 
forehearth changes from 2400°F to 2000°F and across this temperature 
range viscosity varies by a factor of ten. It has been shown that a 

5 

viscosity change of 5 x 10 would affect velocity only by factor of 
three. These facts justify the assumption of constant velocity profiles 
through the channel length. 

To investigate the possible outcome when above assumption is relaxed, 
parabolic profiles are replaced with the solutions of the differential 
equation describing laminar flow in open channels. The differential 
equation to be solved is 



a / au . a / auN dp _ ^ 
ay ay ^ az az^ ' dx ” ‘ 



(3.2) 



Where x is the direction of the flow, u is the velocity, and u is 
the viscosity. Channel width is taken as w and depth being h. A co- 
ordinate system is set at the midstream center with boundary conditions: 

u = 0 @ z = w ( 3 , 3 ) 

u = 0 @ y = h 

^=0 9 z . 0 

^ * 0 9 z = 0 

If viscosity can be considered to be independent of depth and width, 
above equation simplifies into 



2 2 

au.au 1 #_n 



(3.4) 



which is solved with its boundary conditions by the use of infinite series. 
Solution taken from Timoshenko and Goodier [Ref. 7] is 
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(j _ 16Ah 



IT 



2 




n-1 



cosh 2FT 



cosh miw 2h 

W~ 



(3.5) 



Integrating over the area of half-section and dividing by the area (h*w) 
gives the mean velocity F. 



Data for the average mass flow rate gives an estimate for u value, 
which in turn can be used to get a value for a. An experimental value 
of u" from GCIRC report [Ref. 8] was 7.77 ft/hour for green glass. This 
is used in a computer program to evaluate A taking ten terms of the 
series (3.6). Result was 0.03649. With that value, velocity at every 
node of the initial plane is calculated and held constant through the 
channel, as a first try to change previously used profiles. To give an 
idea about the difference obtained in velocity profiles, isovels are 
plotted by 3.0 ft/hour increments in Figure (3b). These changes of main 
computer program increased the running time considerably. Steady-state 
results are reported in Appendix (B), under run (A) for IT = 7.77, run (B) 
for U = 8.54 and run (C) for u" = 7.00. Purpose of the last two runs was 
to get an idea about the effect of mass flow rate on the temperature 
distribution, since the reported value is just a close estimate of the 
actual mass flow rate in the forehearth. Generally better correlation is 
observed to actual data in the center line temperatures, but rather poor 
results observed toward the walls of the channel. Maximum deviation of 
15 degrees F. is calculated in these three runs. 




(3.6) 
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C. VELOCITY PROFILES WITH TEMPERATURE DEPENDENCE 

To build in the viscosity-temperature dependence, available 
viscosity data, Appendix (A), is fitted to a functional form as, 

a+bT 

y = e 

This functional form for viscosity is used to get the values of y at the 
nodes for the numeric calculation and at the same time employing an 
estimated value for (dp/dx). Isovels depending on initial temperature 
distribution are shown in Figure (3c). 

To solve the partial differential equations i.e., heat equation and 
flow equation, the following route is taken: Start with initial tempera- 

ture distribution which is obtained by the use of linear interpolation 
of corner temperatures and calculate velocities at the nodes. Based on 
these velocities, temperatures in the next plane in flow direction are 
calculated by the use of the computer program already developed. All 
the physical parameters are taken as they were in the previous parabolic 
constant velocity profile computations. Computed temperatures at the 
nodes are used to recompute the velocity distribution, which is compared 
with the first velocity array until they agree within the small value e. 
Iterative technique is used in the following way: the first velocity 

array is stored in BA (i,j) array for comparison purposes; after the 
computation of temperatures in one step down the flow direction, they 
are used to get the new velocity array BIJ (i,j). If those two do not 
agree when compared at every node of the plane, a linear combination of 
the two arrays is taken as / 

BIJ(i,j) = W^-BA(i,j) + (1 - W^)BIJ(i,j) ■ 
where W is a weight factor. The BIJ(i,j) array is stored in BA(i,j) 

a 
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Figure 3. Isovels for Channel Flow 
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as the new velocity array, and above cycle is repeated until convergence 
is obtained. Flow diagram is shown in Appendix (D). Numerous runs for 
the first two planes proved that the above method is appropriate. Value 
of W is varied to have an idea of its influence on the convergence and 

di 

also on time to run. Final decision was W, = 0.25 since it cuts the 

a 

running time by half when compared to W = 0.75. 

a 

Due to the inevitable increase in the computer time with the above 
scheme, instead of updating the velocity array at every A x increment 
this is done at every 50 A x increments. This is updating the velocity 
array ten times through the channel length. By the trial runs to assign 
a value to e, it is found that it takes more than 18 iterations for the 
convergence with e = 0.001. Five iterations are needed when e = 0.01 
and 216 seconds to calculate the temperatures at x = Ax and x = 2 ax. 

Taking these facts into consideration and also considering the reliability 
of the data, choice of e = 0.1 is considered appropriate. For run (D) 
velocity array is corrected with 50 A x increments, e value was 0.1, 
running time 1/4 hour. Results at the nodes where data was taken was 
comparable to the best run obtained with the constant velocity calculations 
but not any better. In subsequent runs (D) and (E) the estimate of (dp/dx) 
is increased by 10% and 15% respectively, giving closer fit to the actual 
data. The absolute average deviation was around 4.6 degrees F. in all 
above runs. The largest deviation observed was at y * 5.88 inches and 
z = 8.63 inches with values around -15 degrees F. This deviation was 
also present in the results of the previous model. Another variation 
tried to compute velocity array using 25 A x increments. Results are 
reported under run (F). No improvement is observed but running time is 
increased to twenty minutes. 
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In run (G), boundary temperature at x = XL, y = d and z = w was 
decreased from 2078°F to 2035°F, and velocity updated with 50 A x 
increments. Thus run gave the smallest absolute average deviation 
4.2°F between data point temperatures and predicted temperatures. All 
above runs compute velocity array at the nodes using 10 terms of the 
defining series equation (3.5). When 5 terms of the series is employed 
running time decreased more than three minutes. The last run in these 
lines was run (H) with 5 terms for velocity calculation, e value of 0.1, 
and radiating boundary temperature at the channel corner is increased to 
2140°F from originally used 2125°F. It took 11.8 minutes to run, maximum 
deviation was 12.1°F being 1.8°F less than the best result obtained with 
the previous model. Absolute average deviation from data was 4.22°F. In 
run (I), previously used view factor is replaced with the new one, 
accounting for the reradiation from the channel walls. Maximum deviation 
from the data was 10.6°F, and absolute average deviation was 4.7°F. The 
predicted temperatures at data points for all above runs are listed in 
Appendix (B). 
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VI. PHYSICAL PROPERTY-PARAMETER ESTIMATION 



Physical properties used in the equations must be estimated or 
calculated. Many of these properties and parameters are set at values 



used in the previous model. Since 


this work 


is based 


on the green gl 


only, important parameters that were used in 
kind of glass are listed below. 


original 


model for that 


Radiation Conductivity 


k' 


10.0 


BTU/HR.FT. F 


Density 


p 


147.77 


LBS/FT. 


Specific Heat 




0.375 


BTU/LB. F 


Radiant Surface Emissivity 


^s 


0.9 




Glass Surface Emissivity 


^G 


0.9 




Convective Heat Coefficient 


h 


20.0 


BTU/HR.FT. F 



In addition to the above, parameters a and b in the functional form 
y = are used to get viscosity at the nodes as a function of tem- 

perature. Calculated values for different kinds of glass are given below. 
Plotted viscosity-temperature curve for green glass is presented in Figure 
(4). 

Green 9.07 2.82 10'^ 

Flint 8.60 2.60 10"^ 

To estimate (dp/dx) value; first A is calculated from the average velocity 
data based on the mass flow rate, then this value is used together with 
the average viscosity, since A is defined as 

A = (dp/dx) 

Calculated value for (dp/dx) is 6.849 10"^. 
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Figure 4. Viscosity-Temp Curve. 
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In view factor calculations, which account for the reradiation 
from the channel walls, areas of the radiating surfaces are assumed to 
be equal. Emissivity values are not changed. Geometric view factor 
is taken from Chapman [Ref. 3] as 0.56. 
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V. DISCUSSION 



Comparison of the results obtained by the previous and modified 
models made it clear that the assumption of constant velocity distri- 
bution through the channel does not have significant effect on the 
final results. Predicted temperatures at the nodes were still off by 
small values. Model generally proved flexible and stable for all 
variations made of velocities and of boundary conditions. It still 
does not account for the possible viscosity changes due to pressure 
changes, since pressure changes with position. However, this is not 
considered significant since putting in viscosity-temperature depen- 
dence proved that initial velocities at the nodes decrease by 20-60% 
depending on the position. Further work could involve putting in 
functional forms for density, heat-capacity and thermal -conductivity 
to build in the temperature dependence of these parameters. If this 
is done, the model will be further complicated and limit its possible 
use in a closed loop computer control system. To serve for control 
purposes, time consumption should be decreased considerably. This will 
be possible only by the use of some simplifications. If the accuracy 
in positioning of thermocouples and reliability of the temperature data 
is considered, the present form of the model gives rather good agree- 
ment for the steady-state temperature data. It can possibly be used for 
design purposes. Requirement for a powerful computer is also a draw- 
back for control purposes. An alternate route might be use of a hybrid 
computer instead of a digital ope to make use of the fast integration 
abilities of analog computers to decrease the running time. 
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APPENDIX (A) 



GREEN 


GLASS VISCOSITY DATA 


Temperature°F 


Viscosity (Centi poises) 


2330.4 


316.2 


2292.9 


398.1 


2255.5 


501.2 


2221 . 1 


631.0 


2185.7 


794.3 


2153.1 


1000.0 


2121.2 


1259.0 


2089.4 


1585.0 
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Comparison of data and results predicted with models for green glass 
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Thermocouple Position (in.) 
relative to (0,0,0) point 
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APPENDIX (C) 



FORTRAN VARIABLE 


DEFINITION 


N 

I 

J 

XL 

YL 

ZL 

XEND 

RHO 

CP 

TKP 


Maximum value of X direction index n 

Maximum value of Y direction index i 

Maximum value of Z direction index j 

Length of channel in X direction 
Length of channel in Y direction 
Length of channel in Z direction 
Value of X at which to terminate integration 
Density of the glass 
Specific heat of the glass 
Glass equivalent thermal conductivity in 
Y direction 


TKPP 


Glass equivalent thermal conductivity in 
Z direction 


SIGMA 

H 

ES 

EG 

TSI 

TOGO 

TOOW 

TODO 

TODW 

TLDW 

TLOW 

TA 

TC 


Stefan-Boltzmann Constant 

Convective heat-transfer coefficient 

Emissivity of ceramic radiating surface 

Emissivity of the glass 

Temperature of radiating boundary @ X = 0 

Initial plane corner temperatures 

Corner temperatures @ x = XL 

Temperature of the ambient gas 
Temperature of radiating surface along 
channel centerline 


VMAX 


Max. glass velocity for parabolic 
distribution 


TODWR 


Radiating surface Tern. O x = 0, y = YL, 
z = ZL 


TLDWR 


Radiating surface Tern. 0 x = XL, y = YL, 
z = ZL 


IRNT 


Counter for steady-state program printout 
control 


NTERM 


Number of terms to used in defining series 
for velocity calculation 


ITR 

ITRA 


Counter to recalculate the velocity array 
Counter for number of iterations for 


ITRB 


velocity convergence 

Counter used for non-converging velocities 
at the nodes 


NUP 

WA 

T(i ,j) 


Counter to update the velocity array 
Weight factor for linear combination 
2-Dimensional temperature array 
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BIJ (ij) 


2-Dimensional array for glass velocity 
(V ) in X direction at (i,j) points 
in^y.z plane 


BA (kj) 


2-Dimensional array for intermediate 
velocities 


TAS (ij) 


2-Dimensional array to store tern. @ the 
former plane 
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APPENDIX (D) 

FLOW GRAPH TO CALCULATE VELOCITIES 

ITR = 0 (Set counters to zero) 
ITRA = 0 
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INPLT-OUTPUT SECTION 
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WRITE (6, 91 TAUSBIKI ,TSBSS(K ) , C HGSB ( K J , DL Y SB ( K ) 

3 CONTINUE 

INITIALIZATION AND CALCULATION OF VARIOUS SYSTEM CONSTANTS AND PARAMS 



X 

u 

UJ 

o 



o 







►- 


Q-rvj 
















aoo 














UJ 


U'v 














XQ 


















0X0 














ujin 


to: •• 














o • 




oa 












'vO 




• X 














>--JUJ 


•Hf- 






Q 






^X 


OliJO 


1 V. 






it 






0J< 


UJO 


o> 






O 






OX 


Q>v> 


UJ«J 






* 






» > 




'vUJ 






fO 






X* 


>,JLiJ 


oo 






> 






<^-J-iUiO 


• # 






o 






X>KiujQ •» 


o 


0. 




1 












X 


o 


o 






# -J- 


J'VH-J 




K 


# 


t 








LULLH 


'v 


o 


•H 


h- 




>'v'V-|UJQ 




X 


« 




-J 




'v^fsJUJQ 


(\jO<< 


* 


OJ 


if 


LU 


X > 




O •XK 


>- 






o 


-J ^ 


•J.JUJUJ* tnro 


rsJ(M 1 ^O-lfr 




•> U-5 


H 



■N, LU UJ uJujQQin uu 

Q Q Q QO#* *0 O 

2 'V 'v >ps40* 'C** # > I 

^ ^ ^ «0^<0 

• I I 4> 4 2K2Xi-h>-)ISJ UJUJO^Qu • • • •0«H «UJi--l • 
II II II II II II II QQXXujCMCMrgcg • u 
II II II II II ll»-KXX»Nlfsl II II IIHI-I-II II II lli-IUII II II II 
j II II II i-^<<cDco 

X2£Za.Q.UJUJUJUJUJUJUJLb>>>*-i(\tQ^^(M<^CNJ<^QC<^CM^<Ni 
i^O«-^"5^->QOQQQQQOUUUUUXOUUO<ULUUOO 







X ^ > 


-> 


II 


0«H II 




X 


# M ^UJUJ 


CO 




•HO OOO 




fO 


> --zz 


UJ 


ro 




1- 


kfi 


II 1 


< 


j-jzz 


o 


O II HO II •-•OO 


o 


OO>OOC0oO 


o 










CMro 


< 




corn 


o 




mm 


o 







U 



52 



CALCULATE DENI AND DENJ 
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